The Distribution of Pressures in a Supernova-Driven Interstellar 

Medium. I. Magnetized Medium 

Mordecai-Mark Mac Low 1 
Dinshaw S. Balsara 2 

Jongsoo Kim 2 ' 3 
Miguel A. Avillez 1 ' 4 

ABSTRACT 

Observations have suggested substantial departures from pressure equilibrium 
in the interstellar medium (ISM) in the plane of the Galaxy, even on scales under 
50 pc. Nevertheless, multi-phase models of the ISM assume at least locally 
isobaric gas. The pressure then determines the density reached by gas cooling 
to stable thermal equilibrium. We use numerical models of the magnetized ISM 
to examine the consequences of supernova driving for interstellar pressures. 
In this paper we examine a (200 pc) 3 periodic domain threaded by magnetic 
fields. Individual parcels of gas at different pressures reach widely varying 
points on the thermal equilibrium curve: no unique set of phases is found, but 
rather a dynamically-determined continuum of densities and temperatures. A 
substantial fraction of the gas remains entirely out of thermal equilibrium. Our 
results appear consistent with observations of interstellar pressures. They also 
suggest that the high pressures observed in molecular clouds may be due to 
ram pressures in addition to gravitational forces. Much of the gas in our model 
lies far from equipartition between thermal and magnetic pressures, with ratios 
ranging from 0.1 to 10 4 and ratios of uniform to fluctuating magnetic field of 
0.5-1. Our models show broad pressure probability distribution functions with 
log-normal functional forms produced by both shocks and rarefaction waves, 
rather than power-law distributions produced by isolated supernova remnants. 
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The width of the distribution can be described quantitatively by a formula 
derived from the work of Padoan, Nordlund, & Jones (1997). 

Subject headings: Turbulence, ISM:Kinematics and Dynamics, ISM:Magnetic 
Fields 
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1. Introduction 

Theoretical models of the interstellar medium (ISM) have generally followed Spitzer 
(1956) in assuming that the interstellar gas is in pressure equilibrium. Field, Goldsmith 
& Habing (1969; hereafter FGH) demonstrated that the form of the interstellar cooling 
and heating curves for temperatures below about 10 4 K allowed a range of pressures in 
which two isobaric phases could exist in stable thermal equilibrium. Although the details 
of this model turned out to be incorrect, due to their assumption of a cosmic ray flux much 
higher than subsequently observed, Wolfire et al. (1995) demonstrated that photoelectric 
heating of polycyclic aromatic hydrocarbons or other very small dust grains could recover 
a two-phase model for the neutral ISM. McKee & Ostriker (1977; hereafter M077) used a 
similar framework for the cold ISM, but incorporated the suggestion by Cox & Smith (1974) 
that supernovae (SNe) would create large regions of hot gas with T ~ 10 6 K, to produce 
a model of a three-phase medium. They were the first to relax the assumption of global 
pressure equilibrium, noting that SN remnants (SNRs) would have varying pressures. They 
assumed only local pressure equilibrium at the surfaces of clouds in order to determine 
conditions there. 

Measurements of the ISM pressure were performed using Copernicus observations of 
ultraviolet (UV) absorption lines from excited states of C I by Jenkins & Shaya (1979) 
and Jenkins, Jura & Loewenstein (1983). They found greater than an order of magnitude 
variation in pressures in the cold gas traced by the C I, with most of the gas having 
pressures P/k < 10 4 K cm -3 , but a small fraction reaching pressures P/k > 10 5 K cm~ 3 . 
These findings have recently been confirmed and extended using Space Telescope Imaging 
Spectrograph (STIS) observations by Jenkins & Tripp (2001). Bowyer et al. (1995) & 
Berghofer et al. (1998) compared average pressures derived from observations with the 
Extreme UV Explorer of shadows cast by clouds of neutral hydrogen within the Local 
Bubble with the pressure of the warm cloud surrounding the Sun and found a difference of 
as much as a factor of 25. McKee (1996) argued that these estimates were too extreme, but 
nevertheless concluded that the observations showed a pressure variation of at least a factor 
of five. 

In this series of papers we show that dynamical models of a SN-driven interstellar 
medium do not produce an isobaric medium controlled by thermal instability, but rather a 
medium with a broad pressure distribution controlled by the dynamics of turbulence, with 
rarefaction waves being as important as shocks in setting local pressures. Although the 
cooling curve may determine local behavior, the pressure of each parcel is set dynamically 
by flows primarily driven by distant SNe. The theory of the probability distribution 
function (PDF) of compressible turbulence described by Padoan, Nordlund & Jones (1997; 
hereafter PNJ97) seems able to describe the pressure PDFs we find here. 

Gazol et al. (2001) suggest that turbulence driven by ionization heating results 
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in nearly half of the gas lying in thermally unstable regions, in agreement with the 
observations of Heiles (2001). Earlier dynamical models of the ISM have been reviewed 
by Mac Low (2000). Models including SN driving were first done in two dimensions 
by Rosen & Bregman (1995), and have since been done by a number of other groups. 
Vazquez-Semadeni, Passot & Pouquet (1995) briefly discussed the pressure structure in a 
two-dimensional hydrodynamical model including only photoionization heating. In three 
dimensions, Avillez (2000) used an adaptive mesh refinement code to compute the evolution 
of a 1 x 1 x 20 kpc vertical section of a galactic disk, while Korpi et al. (1999) used a 
single-grid magnetohydrodynamical (MHD) code to model an SN-driven galactic dynamo. 
Slyz et al. (2005) used a single-grid hydrodynamical code to compute supernovae in a 
periodic box including self-gravity and a star-formation recipe. Both of these papers briefly 
discuss pressure structure as well. All of these computations show that the interactions 
between SNRs drive turbulent flows throughout the ISM, and that radiative cooling of 
compressed regions produces cold, dense clouds with relatively short lifetimes. 

We here study in detail the pressure distribution in MHD models of a region in 
the plane of the SN-driven ISM, while in a subsequent paper we will study the pressure 
distribution in a stratified ISM. We use the Riemann framework of Balsara (2000) in 
single-grid mode on a 200 pc cube with periodic boundary conditions on all sides, including 
both radiative cooling and heating, as well as isolated SNs, though neither external or 
self-gravity. 

In § 2 we describe the numerical model. We then give three different analytic 
approaches to the question of the pressure distribution in the SN-driven ISM in § 3, derived 
from the work of M077, PNJ97, and Passot & Vazquez-Semadeni (1998; hereafter PV98). 
In § 4 we use our numerical results to demonstrate that the SN-driven ISM is far from 
isobaric, with order of magnitude pressure variations, but that the distribution of pressures 
takes on a log-Gaussian form that can be well described. These results are compared to 
observations in § 5, and the paper is summarized in § 6. 



2. Models 

The calculations presented here were done using the Riemann framework for 
computational astrophysics, which is based on higher-order Godunov schemes for MHD 
(Roe and Balsara 1996; Balsara 1998a,b; Balsara 2004), and incorporates schemes for 
pressure positivity (Balsara & Spicer 1999), and divergence-free magnetic fields (Balsara 
2004). Balsara & Kim (2004) calibrated the method for the problem computed here. 

We solve the ideal MHD equations including both radiative cooling and uniform 
heating, but neither self-gravity nor external gravity. Our computational domain is a (200 
pc) 3 periodic cube, large enough to contain multiple SN remnants, but small enough to 



- 5 - 



maintain some contact with turbulence in real, stratified galactic disks. As we do not 
include stratification, we must view our models as physics experiments relevant to the ISM, 
rather than truly self-consistent models. We use grid resolutions of 64 3 and 128 3 cells. The 
simulations are started with a uniform density of 2.3 x 10 -24 g cm -3 , threaded by a uniform 
magnetic field in the x-direction with strength either 2 fiG or 5.8 /xG, to cover the range of 
fields suggested by observations of the Milky Way disk (Rand & Kulkarni 1989, Beck 2001). 
Behavior that is common to both sets of models can be deduced to be fairly independent of 
the magnetic field. 

For the cooling, we use a tabulated version of the radiative cooling curve shown in 
Figure 1 of MacDonald and Bailey (1981). This curve follows the equilibrium ionization 
cooling curve of Raymond, Cox & Smith (1976) for 10 6 K < T < 10 8 K; the nonequilibrium 
curve for isochoric cooling derived by Shapiro and Moore (1976) from 10 4 K < T < 10 6 K; 
and then is smoothly extrapolated to zero at 10 2 K. It does not incorporate a sharp cutoff at 
10 4 K due to the turnoff of Lja cooling, nor does it include the expected region of thermal 
instability below 10 4 K. We also include a diffuse heating term to represent processes such 
as photoelectric heating by starlight, which we set constant in both space and time. We 
set the heating level such that the initial equilibrium temperature determined by heating 
and cooling balance is 3000 K. Since the cooling time is usually shorter than the dynamical 
time, we adopt implicit time integration for the cooling and heating terms. 

We set up SN explosions by adding 10 51 erg of thermal energy into a sphere with 
radius 5 pc. The SNe explode at random positions. Unlike Korpi et al. (1999), we allow 
SNe to explode in regions of any density to account for the dominant population of B stars 
(McCray & Kafatos 1987) that last far longer than their parent molecular clouds (Fukui et 
al. 1999). This allows SNe to explode within the shells of pre-existing remnants. In our 
model they are not focussed into associations, however. Some properties of these models 
were described in Balsara, Kim, & Mac Low (2001). 

To test the properties of the code, we set up an isolated SN explosion in an 
unmagnetized medium with the properties of the test problem described by Cioffi, McKee, 
& Bertschinger (1988). We ran this model using the current code and cooling curve until 
the shell hit the edge of our 200 pc box. Figure 1 shows tests with 64 3 , 128 3 and 256 3 zones 
demonstrating that already with 3.13 pc resolution, an isolated explosion agrees excellently 
with the analytic description of Cioffi et al. (1988). 

Table 1 gives the parameters of our runs, including SN rate in terms of the Galactic 
rate, grid resolution Ax, and initial uniform field strength Bq. We take the Galactic SN 
rate, scaled to the area of our grid, to be (1.2 Myr) -1 . The final rms magnetic field < B >, 
and the rms strength of the disordered field (SB) are listed along with their ratio, which is 
discussed in § 4.2. We also list rms Mach number and variance of pressure in our runs for 
comparison with analytic descriptions of turbulence in § 4.3. 
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The evolution of the system is determined by the energy input from SN explosions and 
diffuse heating and the energy lost by radiative cooling. We follow the simulations until 
the total energy of the system, as well as the energy in the thermal, kinetic and magnetic 
energies, has reached a quasi-stationary value for several million years. Similar runs with 
much weaker initial magnetic fields have been used to study the amplification of magnetic 
fields in interstellar turbulence by Balsara et al. (2004). The fields in our runs are amplified 
only modestly, as can be seen by comparing B to (B) in Table 1, because we started our 
runs with relatively strong initial fields. 



3. Analytic Theories 

3.1. Non-interacting Supernova Remnants 

The two-phase theory of FGH assumed pressure equilibrium throughout the ISM, 
with densities and temperatures fully regulated by heating and cooling processes acting 
on timescales shorter than the dynamical timescale of the gas. The introduction of SN 
explosions by Cox & Smith (1974) and M077 required relaxation of the assumption of 
pressure equilibrium, at least inside of expanding SNRs. In particular, Jenkins et al. 
(1983) emphasized that M077 implicitly makes a prediction of the spectrum of pressure 
fluctuations expected from the passage of SNRs expanding in a clumpy medium. 

The theory of pressure fluctuations begins from the scaling in M077, Appendix B, of 
the pressure P in the hot ionized medium, which is given as a function of the probability 
Q(R) of a point being within a SNR with radius at most R, 

I (Q/Q c )- 9 / 14 for Q<Q C 
C \ 0.5(Q/Q c )-°- 9 for Q>Q C , {) 

where Q c is the probability of a point being within a SN remnant large enough for the 
swept-up shell of intercloud medium to have cooled, and P c is the corresponding pressure 
of such a remnant. Note that the immediate loss of pressure after shell cooling produces 
an absence of pressures with 0.5P C < P < P c in this model. For typical values in the 
Milky Way, including a SN rate S = 10~ 13 pc~ 3 yr _1 , M077 find by balancing a number 
of observational considerations that likely values for the parameters in equation (1) are 
Q c = 1/2, and P c /k = 10 3 - 67 cnr 3 K (see M077, eq. 9). 

Equation (1) can then be inverted to give the probability Q(P) of a point being within 
a SN remnant with pressure at least P, as written in Jenkins et al. (1983). The probability 
Q{P) represents a cumulative distribution. The corresponding differential PDF is given by 
—dQ(P). In practice, we compute PDFs by constructing a histogram of pressure values, 
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with bins of finite size AP, so the PDF predicted by M077 is 

^(£)- 23/9 AP for P>P C , 
9P C {kY 19/9aP forP<P c . 



-(dQ/dP)AP={ % V^_ 19/9 _ _ _ (2) 



20Q 



This PDF would diverge towards low pressures if it were not limited by the consideration 
that the pressure in an isolated SN remnant will not fall below the ambient pressure Po, so 
there must be a lower cutoff at that value. (M077 give P /k = 10 310 cm -3 K for the same 
parameters given above.) This model therefore predicts no pressures below the ambient 
value P . 



3.2. Turbulence 

An alternative approach to predicting the distribution of interstellar pressure 
fluctuations can be derived from work on properties of highly compressible turbulence 
by PNJ97 and PV98 5 inspired by the theory of Vazquez- Semadeni (1994). These papers 
considered the PDF of density V(p) in a turbulent, isothermal gas with pressure P = pc? s . 
(PV98 and Nordlund & Padoan [1999] have also considered polytropic equations of state, 
P = Kp 1 where 7 is the polytropic index, but Nordlund & Padoan [1999] note that the 
distribution shifts by less than a factor of two even for polytropic index 7 = 0.3, so we only 
consider the isothermal case here.) By assuming that successive passages of shocks and 
rarefaction waves acting as a random multiplicative process build up density fluctuations 
(Vazquez-Semadeni 1994), they showed that the density distribution is a log-normal 



V{s)ds = — j= exp 
ay2n 



(s - so) 2 



ds, (3) 



2a 2 

where the variable s = In p/po, and p is the mean density of the region. 

The actual value of the variance of the logs of the densities a was extracted from 
numerical models assuming an isothermal equation of state by both groups. Using 3D 
models that have not been completely described in the literature, PNJ97 found 

a 2 = ln(l + M r 2 ms /4), (4) 



5 PV98 suffers from a number of typographical errors produced by a last-minute switch of notation (Passot, 
priv. comm.) to distinguish their scaling parameter M from their rms Mach number M, which we denote 
as M rms . We here enumerate the errors that we are aware of: (1) In the paragraph above their eq. (17), 
M should be used on every occasion. (2) In the text immediately below their Figure 3, a s = ikf, not M. 
(Note however, that their Figure 3 is indeed labelled with M, not M.) (3) There is an extra M in the first 
expression of the middle line of their equation 18, which should be simply u rms /c(s). (4) The first term in the 
bracketed exponential in their equation (20) should contain an additional factor of u^ ms in the denominator. 
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while using very high resolution ID models, PV98 found 

a = M rms , (5) 

where the scaled root mean square (rms) Mach number M rms = v Tias /c s , the ratio of the rms 
velocity to the isothermal sound speed. 

From this formalism, we can derive the equivalent PDF in pressure V(P) by using the 
isothermal equation of state P = pc 2 s . For convenience, we define x = InP, so that 

s — x — In Cg, (6) 

Then the isothermal pressure PDF is 

VWdx = j^j- 2 exp {-^^-} dx, (7) 

where xo = In P(po)- The dispersion of the normalized natural logs of the pressures is thus 
predicted to be the same as that for the pressures, as given in equations 4 and 5. We 
compare our numerical results to these predictions in § 4.3. 



4. Numerical Results 
4.1. Morphology 

We now examine the pressure distribution in our numerical simulations of a SN-driven, 
magnetized, interstellar medium. In Figure 2 we show density, pressure, and temperature, 
as well as magnetic pressure, on cuts through model M2 parallel to the magnetic field. (The 
perpendicular direction appears identical, because the flows are strongly super-Alfvenic, 
with rms Alfven number exceeding four.) Examination of the pressure images immediately 
shows a broad variation in pressures among different regions in all the models, including in 
regions not closely associated with young SNRs. Regions with pressures markedly lower 
than ambient are apparent. 

High temperature regions lie inside young SNRs, while low temperature regions have 
no uniform density and temperature correlation. The highest density regions have sizes of 
dozens of parsecs, and average densities approaching 100 cm -3 , typical of giant molecular 
clouds. This is consistent with the suggestion by Ballesteros-Paredes et al. (1999a) and 
Ballesteros-Paredes, Vazquez-Semadeni & Scalo (1999) that at least smaller molecular 
clouds are formed and destroyed by the action of the interstellar turbulence. Our models do 
not, however include self-gravity as do Li, Mac Low, & Klessen (2005) or Slyz et al. (2005), 
though neither includes magnetic fields (for recent reviews see Elmegreen 2002; Mac Low 
& Klessen 2004), isobaric thermal instabilities (Hennebelle & Perault 1999, 2000; Burkert 
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& Lin 2000; Vazquez-Semadeni, Gazol, & Scalo 2000; and Gazol et al. 2001), or thermal 
conduction (Koyama & Inutsuka 2004). 

Even at SN rates four times those characteristic of the Milky Way, the hot medium 
in our models does not surround discrete clumps of cold and warm gas. Rather, discrete 
regions of hot gas are formed, occasionally intersect, and then seem to be dynamically 
mixed back into the warm gas that fills a substantial fraction of the space. This large-scale 
turbulent mixing appears to enhance the cooling rate, while sheets and filaments confined 
by nearby SNRs seem more effective at slowing down the expansion of SNRs than the 
isolated, discrete clumps of the same equivalent density suggested by M077. These results 
are consistent with the results of Korpi et al. (1999), Avillez (2000), Avillez & Breitschwerdt 
(2004), and Slyz et al. (2005). 

4.2. Thermodynamic Relations 

The first theories of the multi-phase ISM, such as FGH, postulated an isobaric medium. 
Since then, multi-phase models have commonly been interpreted as being isobaric, although 
M077 and Wolfire et al. (1995) actually assume only local pressure equilibrium, not global, 
and M077 did consider the distribution of pressures, as described above in § 3.1. In 
multi-phase models, the heating and cooling rates of the gas have different dependences on 
the temperature and density, so that the balance between heating and cooling determines 
allowed temperatures and densities for any particular pressure. This balance can be shown 
graphically in a phase diagram, showing, for example, the allowed densities for any pressure 
(FGH; for a modern example, see Fig. 3(a) of Wolfire et al. 1995). 

In Figure 3, the thermal-equilibrium curve for the heating and cooling mechanisms 
included in the models is shown as a black line. Only a single phase is predicted at high 
densities as our cooling curve did not include the physically-expected unstable region at 
temperatures of order 10 3 K (Wolfire et al. 1995). Thus, if our model produced an isobaric 
medium, it would be expected to have a single low-temperature phase at uniform density 
given by the point at which the thermal-equilibrium curve crosses that pressure level. 
(Effectively, our cooling curve allows the hotter two of the three phases proposed by M077.) 

Figure 3 show the actual density and pressure of individual zones in the model. Many 
zones at temperatures in the range of 10 3 -10 4 K do lie on the thermal equilibrium curve, 
but scattered up and down it at many different pressures and densities, with no well-defined 
phase structure. Furthermore, a substantial fraction of the gas has not had time to reach 
thermal equilibrium at all after dynamical compression. It appears that pressures are 
determined dynamically, and the gas then tries to adjust its density and temperature to 
reach thermal equilibrium at that pressure. Most gas will land on the thermal equilibrium 
curve when dynamical times are long compared to heating and cooling times. This still 
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leads to all points within the range of pressures available along the thermal equilibrium line 
being occupied, rather than the appearance of discrete phases. Unstable regions along the 
thermal equilibrium curve (Gazol et al. 2001) and off it will also be populated, as observed 
by Heiles (2001) for colder gas, but not as densely, as gas indeed attempts to heat or cool 
to a stable thermal equilibrium at its current pressure. 

Comparison of the models shown in Figure 3(a)-(d) shows that the amount of gas 
that has not yet reached thermal equilibrium depends strongly on the SN rate, but that 
the basic conclusion that gas ends up in thermal equilibrium at a wide range of pressures 
does not. Magnetic fields play a secondary role, with higher field strength increasing the 
total pressure, but somewhat reducing the range of thermal pressures occupied. A stronger 
magnetic field does produce regions of low thermal pressure balanced by magnetic pressure 
that can be seen on and below the thermal equilibrium curve, as we further discuss below. 

The range of pressures observed in our simulations is broader than the region shown 
by Wolfire et al. (1995) to be subject to thermal instability. Even models that included 
a proper cooling curve would produce some gas at pressures incapable of supporting a 
classical multi-phase structure. The mixture of different pressures will, however, produce 
gas at both high and low densities, as well as a smaller fraction of gas at intermediate 
densities that has not yet reached thermal equilibrium. 

Under what conditions are the dynamical times indeed long compared to heating and 
cooling times? We can compare rough analytical estimates of each. The dynamical time is 

^dyn = L/v rrns , (8) 

where L is a characteristic length scale, while the cooling time 

t cool = E/E = kT/nA, (9) 

where E is the thermal energy, and E = n 2 A(T) is the cooling rate as a function of 
temperature T. 

In model M2, for example, the rms velocity volume-averaged over gas of all temperatures 
is f rms = 55 km s _1 . (This is the model with the smallest value of f rms .) If we take 
typical dynamical length scales of L ~ 10 pc, then t dyn ~ 0.2 Myr. We tabulate t coo i 
from our cooling curve in Table 2, normalized to a density of n — 1 cuT 3 . (We note that 
model cooling times in gas at temperatures of 10 3 K ^ T <^ 10 4 K are substantially less 
than physical values, as the MacDonald & Bailey cooling curve does not drop abruptly 
at 10 4 K when Lja cooling shuts off. However, the consequence of this is merely that 
our model overestimates the amount of gas that has reached thermal equilibrium: the 
scatter plot shown in Figure 3 should be even more uniformly filled.) At all temperatures 
10 3 K< T < 10 7 K, we find that t coo i *C td yn , especially at densities n > 1 cm -3 . This agrees 
with our description of gas dynamics setting the local pressure and thermal equilibrium 
only then determining the density and temperature. 
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The separation between dynamical and thermal timescales also sheds light on the study 
by Vazquez- Semadeni et al. (2000) on the effects of turbulence on thermal instability. That 
study suggested that turbulence erases the effects of thermal instability on the interstellar 
medium. Our models show that turbulence probably didn't erase the thermal instability 
entirely, but that the instability only acted locally, under the conditions set for it by the 
larger-scale turbulent flow. Also, thermal instability becomes less important in determining 
the overall distribution of pressures and temperatures when much of the gas has not even 
reached thermal equilibrium. More gas will lie in thermally stable regions than thermally 
unstable ones, but the wide range of available pressures and the lack of complete thermal 
equilibrium in many regions still results in a wide range of properties, despite the nominal 
action of the thermal instability. 

The lower boundary of the heavily occupied region in the pressure-density plane 
appears to be determined by the polytropic behavior of the gas near this cutoff. Fitting 
to its slope yields values of the polytropic index 7 ~ 0.6-0.7. This can be understood as 
due to a cooling curve increasing as a power law (roughly T 2 ' 9 can be fit) balanced by 
heating. Points lying below and to the right of this cutoff line at low pressures and high 
densities appear to be magnetically supported, dropping to very low thermal pressures at 
intermediate densities. This gas does not have low total pressure, only low thermal pressure, 
as shown in Figure 4(6). 

In Figure 4: (a), the relative strength of thermal and magnetic pressure is shown at one 
time for a magnetized simulation. Equipartition of thermal and magnetic pressure (plasma 
(3 = -P g as/-fmag = 1) falls along the solid line. Gas at all temperatures often falls far from 
the equipartition line. In this model, both cold and warm gas have values of (3 as low as 0.1, 
while at the high end, cold gas reaches beta ~ 100, and warm and hot gas have regions with 
(3 as high as 10 4 . Magnetically supported regions occur in both cold and warm gas, as do 
thermally dominated regions. Hot gas can be seen, on the other hand, to be dominated by 
thermal pressure, with low magnetic pressures. In Figure 4(b), the total pressure is shown 
as a function of density for the same zones. Using total pressure, a clear cutoff at high 
densities and low pressures is now plain to see, showing that the previous scatter of points 
beyond it was due to small magnetically supported regions. 

The wide range of magnetic pressures shown in Figure 4a calls into question the 
assumption of equipartition between magnetic and cosmic ray energy usually used to 
derive magnetic field strengths from radio synchrotron emission (Burn 1966, Sokoloff et al. 
1998, Beck 2001). Large fluctuations in the magnetic field strength lead to overestimates 
of the strength of the regular magnetic field for given strength of polarized synchrotron 
emission. Beck (2001) estimates that fluctuations (8B)/{B) ~ 1 lead to as much as a 40% 
overestimate in field strength or a factor of two overestimate in field energy. To measure 
the size of the fluctuations in our models, we take (B) = B Tms , and the fluctuating field 
(SB) to be the rms value of the field with the initial mean field B subtracted off before 
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taking the rms. The resulting values at the end of each of our runs, along with the ratio 
(SB)/ (B) are reported in Table 1. At lower field strengths, the fluctuations indeed approach 
order unity, and even with the high field strengths suggested by the synchrotron emission, 
(5 B) /(B) > 0.5. These results suggest that synchrotron emission observation may yield 
moderate to severe overestimates of the total field strength. 

Finally, let us directly consider the distribution of pressures in gas at different 
temperatures. In Figure 5 we show scatter plots of thermal pressure against temperature. 
We note the concentration of points at temperatures below 10 4 K, which once again reflects 
gas piling up in thermal equilibrium at a wide variety of pressures. The tilt to the left in 
Figure 5 occurs from the effective adiabatic index 7 < 1, with higher pressure gas typically 
having lower temperatures. Individual SNRs with roughly constant pressures are visible 
as stripes in the higher temperature region of the plot. As they cool, they also expand to 
lower pressures. Higher SN rate produces a broader distribution of temperatures in the low 
field case, where SN bubbles expand to fill more of the total volume. 

4.3. Model Probability Distribution Functions 

We have shown that our models display a broad range of pressures. We can quantify 
this by examining the pressure PDFs, as shown in Figure 6. The first point to note, which 
we will return to below, is that all of our models show roughly log-normal pressure PDFs, 
unlike the power-law distributions predicted by the analytic theory derived by M077 
(see § 3.1). The observed distributions rather more resemble the pressure distributions 
suggested by PNJ97 and PV98 (see § 3.2). Not only is the total distribution broad, even 
at the Galactic SN rate (Fig. 6), but so is the distribution for different components of the 
interstellar medium individually. Furthermore, the typical or median pressure at the center 
of the PDF can vary with temperature, as demonstrated by the PDFs of gas at different 
temperature. 

Our results do not appear to depend strongly on numerical resolution as demonstrated 
in Figure 6(a), although the details of the history of each SNR and the total amount of 
energy radiated away will certainly depend on the resolution, as well as on our neglect of 
the physics of the conductive interfaces between hot and warm gas. The pressure PDFs also 
appear to be stable over time, as shown by the comparison of multiple times in Figure 6(b), 
except for the hot gas, especially at the high-pressure end. Individual young SNRs produce 
discrete bumps that move left towards lower pressures as time passes, eventually merging 
with the overall distribution. 

Figure 6(c) shows that increasing the SN rate produces a broader distribution of 
pressures, as generally predicted by both PNJ97 and PV98, and quantified below. Finally, 
Figure 6(d) shows that increasing the magnetic field strength reduces the amount of hot gas 
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as SN remnants are better confined, and reduces the total amount of gas with low thermal 
pressure. 

Using the results for different SN rates illustrated in Figure 6(c), we can quantitatively 
test how well our results agree with the analytic descriptions of PV98 and PNJ97. We fit a 
simple log-Gaussian commensurate with the form of equation (7). These fits are shown in 
Figure 7 for all of our models. Although the Gaussians fit the peaks of our pressure PDFs 
reasonably well, exponential tails do appear in almost all cases. 

We further include a line corresponding to the power law predicted by M077, as given 
in our equation (2), noting that Figure 7 shows the log of the PDF, reducing the power law 
expected by a factor of unity from 23/9 to 14/9 ~ 3/2. This prediction falls in between the 
values seen for the Galactic rate and the higher SN rate near the peak, and overpredicts the 
high-pressure tail at the Galactic rate, while underpredicting the tail at the higher rate. It 
does not appear to give a general description of our results. We further note that there is no 
hint of a gap in the differential pressure distribution corresponding to the lack of pressures 
just below the critical pressure for shell cooling in M077 (see discussion above below eqn. 

(!)• 

We now compare the value of the dispersion of log pressure a that we find in our 
Gaussian fits to the predictions of PNJ97 and PV98 given in § 3.2. As an approximation to 
x , the log of the average pressure, we take the peaks of the PDFs (the modes). We then 
use the rms velocities of each model, to predict a x . In Table 1 we give these values, and in 
Figure 8 we compare the fits to the widths predicted by the rms Mach number for both 
descriptions. Our results agree reasonably well with the prediction for PNJ97, and disagree 
strongly with that for PV98. 

Ostriker et al. (2001) and Li, Klessen, & Mac Low (2003) have noted the failure of 
PV98 to predict the PDFs for three-dimensional, uniformly-driven turbulence with both 
isothermal and non-isothermal equation of state, suggesting that the difference in results 
between PNJ97 and PV98 is due to the use of one-dimensional runs in PV98. The generally 
good agreement we find here with PNJ97 is perhaps somewhat surprising, given that 
their result was derived for uniformly-driven, isothermal, turbulence, while we include a 
non-isothermal equation of state, supernova explosions, and a magnetic field. All of these 
effects appear to play only a secondary role, however. Li et al. (2005) show that the 
assumption of an isothermal equation of state appears to be sufficient to explain the gross 
star-forming properties of galaxies, so perhaps our result is not entirely far-fetched. 

The Gaussian form of the PDF used by both PNJ97 and PV98 was inspired by the 
theory of Vazquez- Semadeni (1994) who argued that the density fluctuation spectrum in a 
strongly supersonic flow is built up by a random series of compressions and rarefactions. We 
can directly trace compressions and rarefactions by examining the distribution of V • v in 
our models (Fig. 9). This distribution is generally symmetric around the peak at zero, with 



-14- 



compressions on the negative side and rarefactions on the positive side appearing nearly 
equally common. This supports the idea that the density and pressure fluctuation spectra 
are generated by repeated compressions and rarefactions occurring with equal probability. 

The question of how mass is distributed among regions of different pressure becomes 
important when considering questions such as the potential ram-pressure confinement of 
molecular clouds. The mass distribution might be expected to markedly diverge from the 
volume distribution given by the PDFs shown previously, as most of the hot gas resides at 
very low densities. In Figure 10 we show mass-weighted distribution functions from each 
high-resolution model. The cold gas dominates the mass-weighted PDF, but it is found at 
the same wide range of pressures, with roughly the same peak pressure, as was suggested by 
the volume-weighted PDFs shown earlier. In particular, a substantial fraction of the mass 
in the cold gas lies at pressures five to ten times higher than the average pressure, even in 
the absence of self-gravity. 

5. Comparison to Observations 

5.1. Ionized and Atomic Gas 

Evidence for a broad distribution of interstellar pressures had already been noted 
as early as the work on excited levels of Cl by Jenkins & Shaya (1979) and Jenkins et 
al. (1983). The point was sharpened with the direct comparison of nearby pressures by 
Bowyer et al. (1995), who compared the pressure of the interstellar medium impinging on 
the heliosphere with the pressure of the interstellar medium averaged over a 40 pc line of 
sight as measured with the Extreme Ultraviolet Explorer. Jenkins & Tripp (2001) have 
used the STIS to extend the work on Cl with better resolved data, concluding that the 
pressure varies by over an order of magnitude both above and below the average value in a 
small fraction of the gas. Our models naturally explain these variations in the context of a 
SN-driven interstellar medium. 

Jenkins et al. (1983) compared their results to the analytic theory of M077, described 
above in § 3. They found a substantially greater column density of low-pressure material 
than that predicted, and rather less high-pressure material. As the analytic theory could 
not predict pressures much less than average, while our models show a broad range of 
rarefaction waves producing a log-normal pressure distribution around the mean, the 
low-pressure results appear consistent. Jenkins et al. were observing emission from excited 
states of Cl in low-temperature gas. We find that much or most of the high-pressure 
gas resides in the hot medium, so, although we do predict a greater total volume of 
high-pressure gas than was predicted by M077, we expect that most of it would have been 
ionized, and therefore not observable by Jenkins et al. (1983), explaining their low column 
densities of high-pressure material. 
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Bowyer et al. (1995) and Berghofer et al. (1998) derived a lower limit to average 
pressure derived from extreme ultraviolet emission along lines of sight to H I clouds in the 
Local Bubble. The cloud distances were determined using photometry of stars superposed 
on the cloud boundaries. They found a pressure of P/k >13, 500-16, 500 K cm -3 for the 
hot medium, with the error depending mostly on the choice of plasma code used to derive 
the emission. They compared this to the value for the pressure in the Local Cloud of 
700-760 K cm~ 3 derived by Frisch (1994) from scattering of solar He I 584Aradiation from 
helium flowing in from the cloud through the heliosphere. McKee (1996) argues that a 
more careful treatment of the unknown ionization fraction could lead to a local pressure 
a factor of three higher, reducing, but not eliminating, the discrepancy. We find greater 
than order of magnitude variations in our models, with pressures reaching values as low as 
a few hundred K cm~ 3 in isolated regions in both sets of models. Interestingly, the lowest 
pressures in our models occur at moderately low densities of order 0.1 cm -3 , while the 
density derived for the Local Cloud by, for example, Quemerais et al. (1994) is 0.14 cm -3 . 

An example of a region with high and low pressure regions intertwined in the manner 
suggested by the EUV observations is the superbubble seen in the corners of Figure 2, 
particularly the part in the upper left corner (note that we are using periodic boundary 
conditions). In this region, even local pressure equilibrium is lacking on scales of tens of 
parsecs, in contrast to all multi-phase models. Similar regions form regularly over time. 

Jenkins & Tripp (2001) found that their results implied an effective polytropic index in 
the cold gas of 7 > 0.9, somewhat higher than the 7 = 0.72 derived by Wolfire et al. (1995) 
for this gas. The suggestion advanced to explain this is that the regions being compressed 
may be smaller than the cooling length scale, and so may begin to behave adiabatically. An 
alternative explanation may be drawn from the broad range of pressures at which gas cools 
in our model: the (relatively sudden) pressurization may happen prior to cooling, rather 
than to already cooled gas, as required by the derivation of 7 by Jenkins & Tripp. 

5.2. Molecular Gas 

Molecular clouds are observed to have broad linewidths suggesting that they are 
subject to pressures as high as 10 5 K cm -3 . Larson (1981) was one of the first to suggest 
that the effective pressure was due to the self-gravity of the cloud. Since then, several 
authors have found that most of the individual clumps in molecular clouds are not in 
hydrostatic balance between turbulent pressure and self-gravity, but rather are confined by 
an external pressure (Carr 1987; Loren 1989; Bertoldi & McKee 1992). The explanation 
offered for this by Bertoldi & McKee (1992) was that the entire cloud was still subject to 
self-gravity, even though individual clumps were not, and so the effects of self-gravity on 
large scales produced pressures that confined the clumps. 
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In our models, we find pressures in high-density regions of order 10 5 K cm~ 3 in the 
absence of self-gravity, as shown in Figure 3. These pressures are sufficient to confine 
observed clumps without invoking self-gravity, suggesting that observed molecular clouds 
may be primarily pressurized by the ram pressure of the turbulent flows in which they are 
embedded rather than being self-gravitating objects. Simulated observations of turbulent 
flows suggest that mass-linewidth relations thought to indicate that they are in virial 
equilibrium may actually be due to a combination of the properties of the turbulence itself 
in the case of the size-linewidth relation, and the properties of the observations in the 
case of the size-mass relation (Vazquez-Semadeni, Ballesteros-Paredes, & Rodriguez 1997; 
Ballesteros-Paredes, & Mac Low 2002). 

Ram pressure is a double-edged sword, however, that can destroy clouds as easily 
as creating them. This is consistent with the suggestion that they are transient objects 
with lifetimes of under 10 7 yr first made by Larson (1981), and emphasized by Ballesteros- 
Paredes, et al. (1999a), Elmegreen (2000, 2002), and Hartmann et al. (2001). Relying on 
ram pressures rather than self-gravity to confine observed clumps in molecular clouds would 
also be consistent with the results of simulations of hydrodynamical and MHD driven, 
self-gravitating, isothermal turbulence that showed self-gravity only acting on small scales, 
with turbulent flows dominating the large scales (Klessen et al. 1999; Heitsch et al. 2001). 
The same flows that confine and destroy the clumps also drive the turbulence observed 
within them, as the background flow stretches, twists, forms, and destroys dense regions 
contained within it. 

Our prescription of isochoric nonequilibrium cooling for 10 4 K < T < 10 6 K, and more 
especially the absence of the physically expected region of thermal instability at T < 10 4 K 
prevents us from having any degree of confidence in the exact amounts of cold gas produced 
at any given pressure. However, our qualitative results appear quite robust, so we expect 
future work to refine the details rather than substantially change the picture of molecular 
clouds forming in transient, high-pressure, high-density regions produced by a supersonic, 
turbulent flow. 



6. Discussion and Summary 

We have examined the distribution of pressures predicted from three-dimensional 
simulations of a SN-driven, magnetized ISM, neglecting vertical stratification and self- 
gravity, but including a distributed heating function. In all the simulations we have run, we 
find pressures distributed over rather more than an order of magnitude around the mean. 
The pressure PDF can be fit well with a log-normal having width given by the analytic 
description of PNJ97 based on three-dimensional simulations of isothermal turbulence, 
and disagreeing with the one-dimensional results of PV98. The log- normal distribution 
is predicted by the theory of Vazquez-Semadeni (1994) that density fluctuations in a 
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supersonic turbulent flow result from a sequence of random compressions and rarefactions. 
We find that compressions and rarefactions are indeed symmetrically distributed, lending 
further support to this idea. The log normal form differs from the power-law distribution 
predicted by M077. The only isobaric regions found in our models are the interiors of young 
SNRs; these are also the only laminar flow regions in an otherwise turbulent environment. 

Most of the mass is in the turbulent gas, however, as is most of the volume, even in 
models driven with four times the Galactic supernova rate. This result is consistent with 
the results of Korpi et al. (1999), Avillez (2000), and Avillez & Breitschwerdt (2004). 
The last group performed a resolution study to establish that numerical diffusion was not 
causing hot gas to cool unphysically. We will study the filling factors of our models in more 
detail in future work. 

We find a broad range of pressures, and a substantial fraction of associated densities, 
far from the thermal equilibrium values. This limits the predictive usefulness of phase 
diagrams based on thermal equilibrium, although thermal equilibrium at the local pressure 
is still the mildly favored state. Gas pressures appear to be determined dynamically. Each 
individual parcel of gas seeks local thermal equilibrium at the pressure imposed on it by the 
turbulent flow. As suggested by previous authors, the phase diagram is only locally valid. 
Isobaric thermal instability in such an environment will lead to regions of the phase diagram 
being mildly disfavored, but no more. Vazquez- Semadeni et al. (2000) and Gazol et al. 
(2001) used two-dimensional simulations with more physically- motivated cooling curves at 
lower temperatures to conclude that the effects of the thermal instability will be barely 
visible in the overall probability distribution functions (but note that these simulations did 
not resolve cooling regions well, so the results should be taken with caution). 

Our results appear consistent with observations that have repeatedly shown the ISM 
not to be isobaric, including those by Jenkins & Shaya (1979), Jenkins, Jura & Loewenstein 
(1983), Bowyer et al. (1995), Berghdfer et al. (1998), and Jenkins & Tripp (2001). 
Although heating and cooling rates remain important for determining the local density 
and temperature, they do not produce a global multi-phase medium with well-separated 
regions of different temperature gas because of the wide range of pressures present and the 
dynamical processing of the gas. Rather, a more continuous distribution of temperatures 
and densities is always present. 

The inference that most molecular clouds must be gravitationally bound because of 
their high observed confinement pressures are called into question by our results. Regions 
with densities approaching the overall densities of GMCs, and pressures an order of 
magnitude above the average interstellar pressure appear in our simulations even in the 
absence of self-gravity. This supports suggestions that star-forming molecular clouds may 
be transient, turbulently-driven, objects (Ballesteros-Paredes et al. 1999ab, Elmegreen 2000, 
2002, Hartmann et al. 2001). Objects that do gravitationally collapse from large scales as 
described, for example by Kim & Ostriker (2001), Elmegreen (2002), or Li et al. (2005) 
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will also form molecular gas, but will quickly form starburst knots in a burst of violent, 
unimpeded star formation. 

Magnetic pressures are also broadly distributed, and do not correlate well with thermal 
pressures. The ratio of thermal to magnetic pressure j3 takes on values ranging from 0.1 
to as much as 100 in gas with temperatures under 10 4 K, and much higher values in hot 
gas. The ratio of fluctuating to mean field in our models (SB) /(B) = 0.5-1, high enough 
to suggest that field strengths derived from observations of polarized synchrotron emission 
may be overestimated (Beck 2001). 

We thank J. Ballesteros-Paredes, D. Cox, E. B. Jenkins, C. F. McKee, T. Passot and 
E. Vazquez-Semadeni for valuable discussions, clarifications of their work cited here, and 
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grant AST 99-85392, while DSB acknowledges NSF grants CISE 1-5-29014, DMS 02-04640, 
and AST 00-98697. JK acknowledges support from the Astrophysical Research Center for 
the Structure and Evolution of the Cosmos funded by the Korea Science and Engineering 
Foundation. Computations presented here were performed at the National Center for 
Supercomputing Applications, which is supported by the NSF, and on the Linux cluster at 
the Korea Astronomy Observatory. This research has made use of NASA's Astrophysics 
Data System Abstract Service. 
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Fig. 1. — Convergence study of shell radius for a single SN remnant in a cold, unmagnetized 
medium. The analytic expression from Cioffl et al. (1988) (line) is compared to simulations 
in a 200 pc box using grids of 64 3 (circles), 128 3 (triangles), and 256 3 (squares) zones. Good 
agreement is obtained at 128 3 resolution. Stair steps in the numerical curves come from only 
being able to measure the size of the shell to the nearest integral zone. 
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Fig. 2. — Two-dimensional slices through the three-dimensional model M8, parallel to the 
magnetic field at a time of 25 Myr, showing density (upper left), thermal pressure (upper 
right), temperature (lower left), and magnetic pressure (lower right). Color bars indicate the 
scale of each quantity, in units of cm -3 for density n, K for temperature T, and dyne cm -2 
for both pressures. The axes are labeled in units of kpc. 
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Fig. 3. — Scatter plot of pressure vs. density at t = 25 Myr in models (a) M2, with low 
SN rate and high field, (b) M4, with high SN rate and field, (c) M6, with low SN rate and 
field, and (d) M8 with high SN rate and low field. We show a subset of 32 3 points sampled 
at intervals of eight points in each direction, at 25 Myr for models M4 and M8, and 100 
Myr for models M2 and M6. Cool gas with logT < 3.9 is shown in blue, warm gas with 
3.9 < logT < 5.5 in green, and hot gas with logT > 5.5 in red. Temperature is given in 
K. The thermal equilibrium curve for the cooling and heating functions in this simulation is 
overlaid as a black line. Note that a wide variation of pressures is seen, and for each pressure, 
different densities occur. 
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Fig. 4. — Scatter plot of (a) magnetic vs. thermal pressure and (b) total pressure vs. density 
at t = 25 Myr in model M8. We again show a subset of 32 3 points. In (a), many of the 
blue cold gas points have been covered by green warm and red hot gas points. The solid line 
shows equipartition between magnetic and thermal pressure. Much of the gas is far from 
equipartition in all temperature ranges. 
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Fig. 5. — Scatter plot of pressure vs. temperature at t — 25 Myr in models (a) M2, with low 
SN rate and high field, (b) M4, with high SN rate and field, (c) M6, with low SN rate and 
field, and (d) M8 with high SN rate and low field. We again show a subset of 32 3 points, 
and give temperature in K. 
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Fig. 6. — Volume-weighted PDFs of pressure for cool gas with logT < 3.9 (blue), warm gas 
with 3.9 < logT < 5.5 (green), and hot gas with logT > 5.5 (red) (a) at different resolutions 
of 1.56 pc (solid, model M8), and 3.13 pc (dotted, model M7), at 25 Myr, (b) at different 
times of 15 Myr (dashed) and 25 Myr (solid) in the 1.56 pc resolution model M8, (c) at 
different SN rates of Galactic (dotted, model M6) at 100 Myr, and four times Galactic (solid, 
model M8) at 25 Myr, and (d) at different initial magnetic field strengths of 2 /iG (solid, 
model M8) and 5.8 fiG (dotted, model M4) at 25 Myr. 
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Fig. 7. — Volume-weighted PDFs (solid) and best fit Gaussians (dotted)of thermal pressure 
for magnetic field strengths of (a) 5.8fiG and (b) for models with resolutions of 1.56 pc 
(red, purple), and 3.13 pc (green, blue), and SN rates of Galactic (blue, purple), averaged 
over 50 times from 75-100 Myr, and four times Galactic (red, green), averaged over 50 times 
from 18.75-25 Myr. Solid lines show the prediction of M077 (taking 14/9 ~ 3/2). Model 
names are noted on the figure. 
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Fig. 8. — Comparison of predictions of PDF width as a function of volume-weighted rms 
Mach numbers from analytic descriptions of supersonic, uniformly-driven, hydrodynamic, 
isothermal turbulence by PNJ97 (solid) and PV98 (dotted) to our low- resolution (crosses) 
and high- resolution (triangles) results. The description of PNJ97 agrees far better with 
our results. Low-resolution models have widths very similar to their corresponding high- 
resolution models, although their Mach numbers differ more significantly in the high SN rate 
cases. 
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Fig. 9. — Volume- weighted PDF of the divergence of velocity normalized by sound speed 
V • v/c s for low field models with supernova rates of (a) Galactic (model M6) and (b) four 
times Galactic (model M8). Negative values greater trace compressions and shocks, while 
positive values trace rarefaction waves and expanding flows. Note the general symmetry 
between positive and negative values at the peaks of the PDFs, supporting the theory that the 
density and pressure distribution is determined by a random series of shocks and rarefactions. 
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Fig. 10. — Mass-weighted PDF of pressure in low field models with SN rates of Galactic 
(dotted, model M6) at a time of 100 Myr, and four times Galactic (solid, model M8) at 25 
Myr, for all the gas (upper left), cool gas with logT < 3.9 (upper right), warm gas with 
3.9 < logT < 5.5 (lower left), and hot gas with logT > 5.5 (lower right). Most of the mass 
is found in cold gas, with a broad distribution around the peak pressure. 
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Table 1. Model Properties. 



model 


T S N a 




Bo 


(B) 


(SB) 


mi 


M rms b 






















Ml 


1 


3.13 


5.8 


6.258 


3.063 


0.489 


0.947 


0.3666 


M2 


1 


1.56 


5.8 


6.443 


3.470 


0.539 


1.039 


0.4100 


M3 


4 


3.13 


5.8 


6.685 


4.430 


0.663 


1.442 


0.6266 


M4 


4 


1.56 


5.8 


7.049 


4.778 


0.678 


1.545 


0.6859 


M5 


1 


3.13 


2.0 


2.436 


1.854 


0.761 


0.947 


0.5222 


M6 


1 


1.56 


2.0 


2.589 


2.071 


0.800 


0.941 


0.5575 


M7 


4 


3.13 


2.0 


2.587 


2.610 


1.009 


1.353 


0.7423 


M8 


4 


1.56 


2.0 


3.071 


3.067 


0.999 


1.580 


0.7378 



a SN rate in terms of the Galactic SN rate 

b Volume-weighted rms Mach number 

c Dispersions of log of pressure derived from fits to PDFs 



Table 2. 


Cooling Times 


logio T 


tcoo\(n/l cm" 3 )" 2 


(K) 


(y r ) 


3 


1.7(5) 


4 


6.0(3) 


5 


1.6(3) 


6 


4.4(4) 


7 


1.2(6) 



